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We analyze the nuclear matter correlation properties in terms of the pair correlation function. 
To this aim we systematically compare the results for the variational method in the Lowest Order 
Constrained Variational (LOCV) approximation and for the Bruekner-Hartree-Fock (BHF) scheme. 
A formal link between the Jastrow correlation factor of LOCV and the Defect Function (DF) of BHF 
is established and it is shown under which conditions and approximations the two approaches are 
equivalent. From the numerical comparison it turns out that the two correlation functions are quite 
close, which indicates in particular that the DF is approximately local and momentum independent. 
The Equations of State (EOS) of Nuclear Matter in the two approaches are also compared. It is 
found that once the three-body forces (TBF) are introduced the two EOS are fairly close, while the 
agreement between the correlation functions holds with or without TBF. 
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I. INTRODUCTION 



The structure and properties of nuclear matter is one of the central issues in the development of nuclear many-body 
theory. Nuclear matter is of great relevance for the physics of supernova, neutron stars and heavy ion collisions, for 
the development of density functionals in nuclear structure studies and for the understanding at fundamental level 
of the low energy baryon-baryon interaction. For a review see reference [1]. Different many-body theories [2, 3] have 
been developed to approach this problem. One can mention the variational method [4-11], the Monte-Carlo method in 
its different versions [12-17] and the diagrammatic expansion methods, in particular the Brueckner-Bethe-Goldstone 
hole-line expansion [2] and the Self Consistent Green Function scheme [18-25]. One of the main goal of this effort along 
the years has been the explanation of the saturation point of nuclear matter that can be extracted phcnomcnologically 
through various experimental method, in particular the analysis of the binding energy of nuclei and of the electron 
elastic scattering cross sections [1]. However, besides the saturation point, one of the most important characteristics 
of nuclear matter is its correlation structure. In fact the presence of a hard core in the nucleon-nucleon interaction 
produces a correlation "hole" between two nucleons that can be described by the correlation function. The latter is 
also determined by the intermediate and long-range interaction, typical of the nuclear two-nucleon potential. The 
correlation function is a key quantity to characterize each many-body scheme and to understand the corresponding 
numerical results. In scattering studies, the spectral function of many fermionic system gives the important quantities 
of interests and the short and long range correlation functions are very important factors for calculating the spectral 
functions [26]. The connection of the correlation function and the spectral function is not straightforward, but it has 
been elucidated in ref. [27] in the Brueckner-Hartree-Fock (BHF) framework and in ref. [28] for the Lowest Order 
Constrained Variational (LOCV) method, where it is particularly transparent. Furthermore several phenomena that 
occur in neutron star matter are closely linked to the correlation function, like e.g. dissipation due to shear viscosity and 
neutrino transport. It appears then natural to look for a comparison between the correlation functions from different 
many-body schemes. In this paper we present a detailed comparison between the Bethe-Brueckner-Goldstone (BBG) 
method [2] and the variational method, as developed within the LOCV framework. Both methods have been applied 
systematically to nuclear matter with different two-body interactions. The results for the saturation point and other 
physical parameters, like the compressibility at high density [29, 30], the critical temperature of the liquid-gas phase 
transition [31, 32], are close but not completely in agreement. One of the main goals of this work is to present an 
analysis of the correlation function that could help understanding the reason of the agreements and the discrepancies 
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by the comparison of the corresponding correlation properties. 

II. THE VARIATIONAL METHOD 

The method of Lowest Order Constrained Variational approach is among the microscopic methods that were 
developed to calculate the bulk properties of homogeneous nuclear fluids such as the saturation quantities by using 
the realistic nucleon-nucleon interaction i.e. Reid68 and A— Rcid (the modified Reid potential with inclusion of 
isobar degrees of freedom) [33] . This method was reformulated to include more sophisticated interactions [34] such as 
UVi4,AVis [35] and charge dependent Reid potential (Reid93) [36]. The LOCV method has been also developed for 
calculating the various thermodynamic properties of hot and frozen homogeneous fermionic fluids such as symmetric 
and asymmetric nuclear matter [37], j3— stable matter [38], Helium 3 [39], electron fluid [40] with different realistic 
interactions. Recently the LOCV formalism was developed for covering the relativistic Hamiltonian with a potential 
which has been fitted relativistically to nucleon-nucleon phase shifts [41]. The LOCV calculation is a fully self- 
consistent technique with state dependent correlation functions. There is no free parameter in this method, except 
those included in the interactions. Considering constraint as in the form of normalization condition is another 
advantage of LOCV formalism. This assumption keeps the higher order terms as small as possible and it also assumes 
a particular form for the long range behavior of the correlation functions in order to perform an exact functional 
minimization of the two-body energy with respect to the short range parts of correlation functions. The functional 
minimization procedure represents an enormous computational simplification over the unconstrained methods, where 
the short range behavior of the correlation functions is parametrized, that attempt to go beyond the lowest order 
[42]. To test the convergence of the LOCV method for nuclear matter and helium 3, the calculations were performed 
beyond the lowest order and the three-body cluster energy was evaluated with both the state averaged and state 
dependent correlation functions [43]. The smallness of the normalization (the convergence parameter) and of the 
three body cluster energy indicated that at least up to the twice empirical nuclear matter saturation density, the 
cluster expansion converges reasonably and stopping after two-body cluster terms is a fair approximation. In the 
LOCV method, we use an ideal Fermi gas type wave functions, fa , for the single particle states and we employ the 
variational techniques to find the wave function of interacting system [33]-[36],i.e, 

* = F$, (1) 

where $ is the uncorrclated Fermi system wave function (Slater determinant of plane waves) and the factor F(l, 2.., A) 
is the many-body correlation function, defined as product of two body correlation functions f(i,j) (Jastrow form) 
and assumes that they are operators, 

F = Sl[f(i,j), (2) 

i<j 

where S is a symmetrizing operator. The many-body energy term E[f], which is a functional of the /' s, is calculated 
by constructing a cluster expansion for the expectation value of Hamiltonian H of the system. 

1 < * I H I * > 

E W = A <*\i> + > E °- & 

where Eq is the true ground state energy and A is the particle number. In the lowest order we truncate the above 
series after E 2 i.e. two-body energy. The one body term E\ is independent of the f and is just the familiar Fermi gas 
kinetic energy. The two-body energy term is defined as, 

E2 = X] < ij ' I W I V ' >a ; I i: > >a= \ij > ~\ ji > 

ij 

W = -|^[/(1,2),[V 2 ,/(1,2)]] + /(1,2)V(1,2)/(1,2), (4) 

and the two-body anti-symmetrized matrix element < ij \ W\ij > a are taken with respect to the single-particle 
functions composing fa i.e. plane-waves. By inserting a complete set of two-particle state twice in above equation 
and performing some algebra we can rewrite the two-body term as a functional of correlation functions [33, 34, 36] 
. In this equation V(l,2) is phenomenological nucleon-nucleon potential such as Reid type, UV14, and AVi$. At this 
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stage, we can minimize the two-body energy with respect to the variations of the correlation functions [33, 34, 36] , 
but subject to the normalization constraint [33]-[41] : 

j < ij\h 2 (l,2) - f(l,2)\ij > a = 1. (5) 
The function h(12) is the modified Pauli function, which for the symmetrical nuclear matter take the following form: 

MllJ) _( 1 .| ( *M ) »)-», (9) 

where ji{r\2) is the well known spherical Bessel function of order 1. Note that [x =< > —1] plays the role of a 
smallness parameter in the cluster expansion. The above constraint introduces a Lagrange multiplier through which 
all the correlation functions are coupled. Then we can write sets of uncoupled and coupled Euler-Lagrange differential 
equations with respect to the correlation functions. The constraint is incorporated by solving these Euler-Lagrange 
equations only up to a certain distance where the logarithmic derivative of correlation functions matches those of 
Pauli function and then we set the correlation functions equal to Pauli function . As we pointed out before, there is 
no free parameter in our LOCV formalism i.e. the healing distance is determined directly by the constraint and the 
initial conditions. 



III. THE BBG EXPANSION 

One of the most known and used microscopic many-body approach to the theory of nuclear matter is the Bcthc- 
Brucckner-Goldstone (BBG) expansion [2]. In this scheme the original nucleon-nucleon interaction is systematically 
replaced by the so-called G-matrix, that describes the two-nucleon scattering amplitude inside the medium. A 
modified perturbative expansion is then developed in terms of this effective interaction and the different terms can be 
represented by diagrams. The G-matrix can be defined also for singular interaction, e.g. with a hard core, and it is 
expected to be "smaller" than the original NN interaction. Although all modern realistic NN interactions introduce 
a finite repulsive core, it is however quite large, and therefore in any case a straightforward perturbative expansion 
cannot be applied. As discussed in the presentation of the variational method, the repulsive core is expected to modify 
strongly the ground state wave function whenever the coordinates of two particles approach each other at a separation 
distance smaller than the core radius c. In such a situation the wave function should be sharply decreasing with the 
two particle distance. The "wave function" of two particles in the unperturbed ground state (j>o can be defined as 
(k l7 k 2 < k F ) 

0(n,r 2 ) - (Ui'iir^lir^auM = ^(k 1 +k 2 )-R e ,(k 1 -k 2 ). r /2 ; (7) 

where £i ^ £2 are spin-isospin variables, and R = (ri + r 2 )/2, r = (ri — r 2 ) are the center of mass and relative 
coordinate of the two particles respectively. Therefore the wave function of the relative motion in the s-wave is 
proportional to the spherical Bessel function of order zero j'o(fcr), with k the modulus of the relative momentum 
vector k = (ki — k 2 )/2. The core repulsion is expected to act mainly in the s-wave, since it is short range, and 
therefore this behavior must be strongly modified. In the simple case of k = the free wave function jo(kr) — > 1, and 
schematically one can expect a modification, due to the core, as depicted in Fig. 1. The main effect of the core is 
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FIG. 1: Schematic representation of the expected effect of the core repulsion on the two-body wave function in nuclear matter. 



to "deplete" the wave function close to r = 0, in a region of the order of the core radius c. Of course, the attractive 
part of the interaction will modify this simple picture at r > c. If the core interaction is the strongest one, then 
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the average probability p for two particles to be at distance r < c would be a measure of the overall strength of the 
interaction. If p is small, then one can try to expand the total energy shift AE due to the interaction in power of 
p. The power p n has, in fact, the meaning of probability for n particles to be all at a relative distance less than c. 
In a very rough estimate p is given by the ratio between the volume occupied by the core and the average available 
volume per particle 

with ^-d 3 = p^ 1 . From Eq. (8) one gets p w -^(kpc) 3 , which is small at saturation, fe^ = 1.36 /to -1 , and the 
commonly adopted value for the core is c = 0.4 fm. The parameter remains small up to few times the saturation 
density. 

The terms of the expansion can now be ordered according to the order of the correlations they describe, i.e. the 
power in p they are associated with. It is easy to recognize that this is physically equivalent to grouping the diagrams 
according to the number of hole lines they contain, where n hole lines correspond to n-body correlations. In fact, 
an irreducible diagram with n hole lines describes a process in which n particles are excited from the Fermi sea and 
scatter in some way above the Fermi sea. Equivalently, all the diagrams with n hole lines describe the effect of clusters 
of n particles, and therefore the arrangement of the expansion for increasing number of hole lines is called alternatively 
"hole expansion" or "cluster expansion". For a pedagogical introduction to the BBG expansion see refs. [2, 44], where 
references to more technical reviews can be found. In ref. [44] the connection of BBG and the variational method 
is discussed. The relation between the two approaches turns out to be more transparent if the BBG expansion is 
reformulated in terms of the coupled cluster method (or e s method) [45] . According to this scheme the wave function 
of the ground state is written 

I*) = e S |*>, (9) 

where S is a correlation operator containing a set of n-body terms which produce excitations of n particles from below 
to above the Fermi sea. This method has also a variational character, in the sense that the variation is performed not 
on the ground state wave function but on these correlation terms [45, 46]. Then a set of coupled equations is obtained 
for the n-body correlation functions. The expansion of this set of equations in terms of the order of the correlations 
is equivalent to a re-ordering of the hole-line expansion in the BBG theory [47]. At the 2-body level of approximation 
the method is equivalent to the so-called Brueckner approximation [44, 47] in the BBG hole-line expansion, and the 
operator S reduces to a two-body operator §2 

S2 = {Kk l 2\S2\kik2)a\k l l )a\k l 2 )a(k 2 )a{k 1 ), (10) 

where the k ' s label hole state , i.e. inside the Fermi sphere, and the k' ' s particle states, i.e. outside the Fermi 
sphere. Each quantity k indicates momentum k and spin-isospin quantum numbers. The function S2 is the so called 
"defect function" of the Brueckner scheme. It can be written in term of the G-matrix and it is just the difference 
between the in- medium interacting and non interacting two-body wave functions [2, 44]. The different terms of the 
summation commute with each other and expanding the exponential in Eq. (9) one gets the product of the correlation 
operators over all sets of momenta, 

I*) = n {fe} [l + Y,(Kk'2\S2\k 1 k 2 )a\k' 1 )aHk' 2 )a(k 2 )a(k 1 )]\^) 7 (11) 

where the product is over all disjoint pairs of momenta k\,k2 1 in a given partition of the set of all momenta, in 
agreement with the Brueckner scheme, which is an independent pair approximation. Higher orders in the expansion 
vanish because they include powers of annihilation or creation operators. In the square bracket one can recognize 
the two-body wave function. After Fourier transformation to coordinate representation and assuming the defect 
function to be local and independent of total momentum, this expression acquires the same form as in the variational 
method, where the two-body wave function plays the role of the correlation factors f(i,j). However there are relevant 
differences with the variational method. First of all the BBG expansion is not explicitly variational, although, as 
already mentioned, one can recast the expansion in terms of the e scheme, which can be formulated by means of a 
particular variational procedure [44, 45]. Second, the G-matrix, and therefore the defect function, is in general highly 
non-local, which means that the two-body wave function is dependent also on the initial momenta in the Fermi sea, 
as well as on the total momentum. This would imply a correlation factor in integral form for the variational scheme. 
Furthermore in BBG expansion one introduces a single particle auxiliary potential, in order to increase the degree 



■5 



of convergence of the expansion. This potential is usually called Brucckner potential and it is determined with a 
self-consistent procedure [2]. In the variational method no single particle potential is introduced in the minimization 
procedure. Of course it is hidden in the mean value of the hamiltonian, but it can be calculated only after the optimal 
many-body wave function and energy have been obtained, by adding a tiny fraction of particle to the system [48]. 
Finally in the variational method the correlation function is introduced in the mean value of both the kinetic energy 
and the interaction term. It is a peculiarity of the BBG expansion that the total energy is written as the sum of the 
unperturbed kinetic energy and the correlated interaction energy. The latter includes of course implicitly the effect 
of the correlation on the kinetic energy due to the momentum dependence of the single particle potential and of the 
G-matrix. 

It is one of the main purposes of this work to explore the consequences of these differences on the correlation 
properties of the ground state. In turn, the study provides a detailed view of the nuclear matter correlations. Since 
for both LOCV and the BBG expansion three-body correlations turns out to be only a fraction of MeV around 
saturation density [44, 49, 50], we restrict the comparison to two-body correlations. In any case the two-body 
correlation functions are determined at the BHF level for the BBG expansion and at the two-body Jastrow-like 
factors for the variational method. 



IV. FORMAL AND NUMERICAL COMPARISON 

In order to formulate a meaningful comparison between the two-body microscopic methods, we introduce a mixed 
representation of the correlation functions. In the expansion of Eq. (11) we separate relative and total momenta and 
perform the Fourier transformation on the momenta k', i.e. the final ones above the Fermi sea. One gets in this way 
the correlation function F in coordinate representation, which is dependent on the initial relative momentum and on 
the total momentum 

/ 0f £lcl ' r[Sqq ' + (<z ' |52(p)k) ] = Fb ^^ p ^ ( 12 ) 

where the defect function can be written in term of the G-matrix 

W\S 2 (P)\q) = ^^ W\G{P)\q), (13) 

where Q is the average Pauli operator and e is the average two particles excitation energy, see the Appendix for more 
details. Because of this averaging, the denominator in Eq. (13) can vanish. The integral of Eq. (12) is meant as 
principal value, in agreement with the BHF calculations of the nuclear matter EOS. 

The correlation function can be expanded in partial waves and one can define a correlation function for each two- 
body channel, identified by the quantum numbers ISJT of the relative angular momentum, total spin, total angular 
momentum and total isospin, respectively. As shown in the appendix, the correlation function Fb has to be compared 
with the corresponding correlation function Fy for the variational method 

F V (r,q) = f(r) ■ ji(qr), (14) 

where ji is the spherical Bessel function of order I and f(r) is the correlation function of e.g. Eq. (5). It is essential 
to notice the factorization of the free wave function, characteristic of the variational method. For the Brueckner 
correlation function Fb this property does not hold, which embodies the non-locality of the G-matrix. However it can 
hold approximately, and this can be verified by e.g. the numerical comparison between the two correlation functions. 
Details on the formal comparison between Fy and Fb can be found in the Appendix. 

We consider symmetric nuclear matter around saturation and we take the potential Argonne v\% [35] as the two- 
body nucleon-nucleon interaction. At the Fermi momentum kp = 1.36 fm , corresponding to density 0.17 fm~ 3 , we 
compare in Fig. 2 the correlation functions Fy(r) and Fb(t) at the relative initial momentum q = 0.1 fm -1 and at 
zero total momentum P. In this case the correlation functions are calculated for the 1 Sq channel. In the variational 
method a small hard core of radius R c — 0.1/m is introduced for numerical reasons, which is apparent from the 
figure since the correlation function Fy is zero below the core radius. Both correlation functions feel in any case the 
repulsive, but finite, core of the interaction and they decrease sharply at short distance. They agree closely above 
the small core radius R c . At large distances both correlation functions reach the expected value of 1, but just above 
the repulsive core they exceed 1, due to the attractive part of the NN interaction. In this region they practically 
coincide. A small discrepancy is observed at intermediate distances, where Fy is slightly larger than Fb- To be more 
quantitative, we calculated the mean absolute deviation for r > 0.25 fm. We found a value below 2%, as in all cases 
we are going to consider in the following. 



6 



The comparison for the So channel, but for q = 0.5 fm -1 , is reported in Fig. 3. In this case already at moderate 
distance the two-body wave function F starts to oscillate since it smoothly merges into the free wave function, i.e. the 
Bessel function (of order in this case). The same agreement between Fy and Fb is observed. This result indicates 
that the factorization of Eq.(14) is approximately valid also for the correlation function Fb of the BBG expansion. It 
is also an indication that the defect function is approximately local. 

Notice that in the numerical calculations the correlation functions are multiplied by r 2 and therefore the contribution 
of the small distances is vanishing small. This is illustrated in the same figures, where the correlation functions 
multiplied by r 2 are reported. In this case the very close agreement is apparent. A similar trend is obtained for the 
3 S± channel, Fig. 4. In the channels with higher partial waves the agreement is even better. To put in evidence 
the tiny differences, we have reported in an amplified scale the correlation functions in Fig. 5 for the l Pi channel 
and in Fig. 6 for the 3 Pi channel. Notice the change of scale with respect to the previous figures. In these cases 
the centrifugal barrier suppresses further the two-body wave functions at short and intermediate distance. At larger 
distance, outside the considered range, the correlation functions merges into the proper Bessel function and then they 
obviously coincide. 

We also checked the dependence on the total momentum that is present in the two-body wave function. It turns 
out that this dependence is quite weak, see Fig. 7, which justifies the assumption, intrinsic in the variational method, 
of neglecting such a dependence. Finally we have introduced the three-body forces (TBF) in the calculations, both in 
the LOCV and the BBG schemes. It is well known that TBF are necessary if the phenomenological saturation point 
of nuclear matter has to be reproduced. At the level of two-body correlation approximation, as BHF and LOCV, the 
TBF arc reduced to an effective two-body force by averaging on the position and on spin-isospin of the third particle 
[51]. The averaging involves the two-body correlation itself. In principle the original TBF can be derived within the 
nucleon-meson model of nuclear forces. This procedure turns out to have only a limited success [52, 53] and requires 
in any case the tuning of the parameters (masses and coupling constants) to get a reasonable saturation point. The 
latter can be obtained only with the Bonn B potential [54, 55] as two-body forces [52, 53]. We prefer to follow a more 
pragmatic point of view. We used the Urbana IX model and treated the TBF according to the method adopted in 
ref. [31], where the averaging is performed by using a schematic two-body correlation function. We then tune the 
(two) parameters of the TBF to get a good saturation point for BHF and we use the same values in LOCV. Around 
saturation the contribution of TBF is relatively small in absolute value, about 1-3 MeV, in comparison with the total 
correlation energy that is about -40 MeV at this density. It is slightly repulsive, and as a consequence the two-body 
wave function is further reduced. This is illustrated in Fig. 8. Since the effect is quite small, as expected by the 
relatively weakness of the TBF, Fig. 8b shows a blow up of the small distance region. It looks that the effect of the 
TBF is slightly larger for the BHF method. 

V. DISCUSSION AND CONCLUSION 

We have studied the correlation properties of nuclear matter both in the variational LOCV method and in the BHF 
scheme. In particular we have shown that one can identify the variational (generalized) Jastrow factor Fy(r) with 
the BHF correlation function Fs{r) = 1 + g(r), where g(r) is the so-called Defect Function. Despite the additional 
total and relative momentum dependence of Fb, not present in Fy, and the different method of approximation, it 
turns out that the two correlation functions are quantitatively quite similar. This is true for each two-body channel, 
with or without the inclusion of the three-body forces. To see the possible relation of the small differences between 
the LOCV and BHF correlation functions to other nuclear matter properties, we have computed the nuclear matter 
Equation of State (EOS) in the two theoretical schemes. The results are reported in Fig. 9. The two lower curves, 
labelled 2BF, correspond to the EOS with two-body forces only, while the two upper curves, labelled 2BF + 3BF, 
correspond to the EOS when the (same) three-body forces are also included. One can notice that the two EOS are 
much more similar when the three-body forces are included. This is in line with the similar finding [56] that the EOS' 
s with different NN interactions become much closer when the (same) three-body force is included. 

In the variational method the average kinetic energy is affected directly by correlations. The total correlation 
energy includes a kinetic energy part and a potential part, see Eq. (4). The breakdown of the two contributions as 
a function of density is reported in Table I for the case where TBF are included. For comparison the total potential 
energy of the BHF calculations is also reported. In the BHF scheme the kinetic energy is not explicitly modified [2] , 
and the whole correlation energy is contained in the potential energy coming from the G-matrix contribution. The 
modification of the kinetic energy is embodied in the momentum dependence of the G-matrix and in the self-consistent 
single particle potential, which also affects the total binding indirectly since it determines the entry energy of the 
G-matrix. From the results it looks that the connection of the EOS and the details of the correlation function is not 
so straightforward. This is apparent if we calculate the correlations functions at twice the saturation density. They 
are displayed in Fig. 10a for the case with only two-body forces. The agreement between the two correlation looks 
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TABLE I: Nuclear matter correlation energy per particle in LOCV and in BHF as a function of the density p. The first 
column (K.E.) for LOCV gives the modification of the kinetic energy due to the two-body correlation, the second one (P.E.) 
the potential part and the third one their sum. For comparison the BHF total correlation energy is reported in the last column 
(BHF). The three-body forces are included. 



P(fm- J ) 


LOCV 


BHF 




K.E. 


P.E. 


TOT 




0.10 


11.24 


-40.75 


-29.51 


-29.24 


0.17 


16.77 


-56.09 


-39.32 


-37.97 


0.20 


18.94 


-61.20 


-42.26 


-40.42 


0.30 


25.62 


-71.18 


-45.56 


-43.60 


0.34 


28.20 


-72.34 


-44.14 


-43.10 


0.40 


32.12 


-71.83 


-39.71 


-40.46 


0.50 


36.68 


-64.32 


-27.64 


-31.07 



insensitive to the introduction of the TBF, see Fig. 10b and indeed quantitatively the disagreement, as anticipated 
before, is below 2% for r > 0.25fm, with or without TBF. Despite small variations can be relevant, it looks unlike that 
this deviation can be considered responsible of the fact that the disagreement between BHF and LOCV is reduced by 
several MeV at this density once TBF are introduced. It has to be noticed that in BHF there is no simple way to relate 
the binding energy to the correlation function, which is not directly involved in the BHF expression for the correlation 
energy. The change of the binding is clearly due to the direct effect of the change in the nucleon-nucleon force due 
to the TBF. The only effect of TBF, on the correlation function, see Fig. 10b, seems to be a very small decrease 
at intermediate distance of Fy with respect to Fb- This could suggest that the good agreement of the EOS is the 
result of a redistribution of the attractive and the repulsive contributions to binding. To make easier the qualitative 
estimate of the relevance of the TBF, we have reported in Fig. 11 the comparison of the correlation functions with 
and without TBF at the densities 0.16 fm -3 and 0.32 fm~ 3 , both for LOCV and BHF. At increasing density the effect 
of TBF increases, but the effect looks larger for BHF. Also in this case no systematic trend is observed in relation to 
the corresponding EOS. It has been found in ref. [57] that also the spectral function has a mild dependence on the 
presence of TBF. Beside the EOS, other quantities, like transport coefficients or neutrino and electron scattering cross 
sections, are probably more directly related to the correlation and spectral function [58]. The analysis of this point is 
left to a future work, but in any case no major discrepancy can be expected between BHF and LOCV schemes. 
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Appendix 

In this Appendix we give some details on the formal comparison between the effective correlation factor Fb that is 
present in the ground state wave function of the Brucckncr approximation, within the BBC hole-line expansion, and 
the corresponding correlation factor Fy in the variational LOCV approximation. The unperturbed ground state $ is 
the anti-symmetrized product of N single particle momentum states 

|$) = U ki al\0), (A.l) 

where \0) is the vacuum state and the ki include spin-isospin variables. Then Eq. (11) can be rewritten 

I*) - n {felfe2}J F felife2 |0) 

(A.2) 

Fk lk2 = E k > lk > [Sk'kJk^ + (Kk> 2 \QS 2 \hk 2 )] at(fci)at(^), 

where the summations in F are over all momenta and we have introduced the Pauli operator Q that restricts the 
momenta k[k' 2 outside the Fermi sphere, while the momenta k\k 2 are inside the Fermi sphere. 

It is convenient to introduce the wave function of the correlated ground state |^) by taking the scalar product with 
the anti-symmetrized N-particle coordinate states 



\nr 2 ....r N ) = U^(n)\0) = \{n}) 



(A.3) 
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where V^( r i) is the creation operator of a particle at the position n (including spin-isospin variables). One gets 

*({*}) - ({nm = A {rt} n {klk2} f klk2 (n ir] ) (A.4) 

where the operator A anti-symmetrizes the N coordinates and 

fhk 2 {n,rj) = [ S K*M>* + (k[k 2 \QS 2 \hk 2 )] (nlki)^), (A.5) 

which is the Fourier transform of the defect function. The variables r, and rj are two generic coordinates among the 
N anti-symmetrized ones. Introducing the coordinate representation for the defect function, one gets 

fk^in^j) = (r-i|fci)(rj|fc 2 ) + (r l r j \QS 2 \k 1 k 2 ) 

(A.6) 

= <ri|fci)fo|fc2> + J d^d 3 ^ {nr.lQS^) {r'^irrih). 

We consider the relative coordinate = (rj — rj) and center of mass coordinate Rij = (r^ + rj)/2 and notice that 
the defect function QS^ is diagonal in the total momentum P. If furthermore we assume that the defect function is 
local, one gets 

fkM^rj) = (rij\q)(Rij\P) + Jd%(r lJ \QS 2 (P)\r' ll )(r' l3 \q)(R l] \P) 

(A.7) 

= [1 + gMinAqHRiAP), 

where q is the relative momentum and 

(r tl \QS 2 {P)\r' l3 ) = g{n 3 )5{n 3 - r[ 3 ). (A.8) 

Here the dependence on the total momentum of the defect function has been neglected. This result shows that, under 
the stated assumptions, the correlated wave function can be written as 

*({ri}) = AU klk2 [1 + g( rij )} (nlhUrjlki) (A.9) 

which has the form of the variational wave function, if we identify the factor 1 +g with the correlation function f(r) of 
the variational method. The defect function in the mixed representation F B has then to be compared with F v (r){r\q), 
as discussed in the text. Both F B and F v can be expanded in partial waves and compared channel by channel. 
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FIG. 2: Correlation function in the 1 So channel for the LOCV and BHF approaches. The same correlation functions multiplied 
by r 2 are also shown. The momentum q — O.lfm -1 is the relative momentum of the two correlated particles. 
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FIG. 7: The correlation functions as in Fig. 4 at different total momentum P of the two correlated particles. 
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FIG. 8: In panel (a) the correlation function in the Si channel is reported, with and without three-body forces. Panel (b) is 
the blow up of the plot in panel (a) within a region at small distances. 




FIG. 9: Binding energy per particle as a function of the density p in symmetric nuclear matter for the LOCV and BHF 
approaches. The two lower curves, labelled 2BF correspond to calculations with two-body force only. The two upper curves, 
labelled 2BF + 3BF correspond to calculations with the inclusion of the three-body forces. 
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FIG. 10: Correlation function in the 1 So channel for the LOCV and BHF approaches at the density p — 0.32 fm~ 3 . The 
same correlation functions multiplied by r 2 are also shown. The momentum q = O.lfm -1 is the relative momentum of the two 
correlated particles. Panel (a) : only two-body forces. Panel (b) : also three-body forces are included. 
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FIG. 11: Comparison of the effect of the TBF on the correlation function at two different densities, 0.16 fm~ 3 (full lines) and 
0.32 fm -3 (dashed lines). At each density the lower curves include the TBF. Panel (a) refers to BHF, panel (b) to LOCV. The 
meaning of the other labels is as in previous figures. 



